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Summary 

A research program is in progress to develop strain rate dependent deformation and failure models for the 
analysis of polymer matrix composites subject to high strain rate impact loads. Strain rate dependent inelastic 
constitutive equations have been developed to model the polymer matrix, and have been incorporated into a 
micromechanics approach to analyze polymer matrix composites. The Hashin failure criterion has been 
implemented within the micromechanics results to predict ply failure strengths. The deformation model has 
been implemented within LS-DYNA, a commercially available transient dynamic finite element code. The 
deformation response and ply failure stresses for the representative polymer matrix composite AS4/PEEK have 
been predicted for a variety of fiber orientations and strain rates. The predicted results compare favorably to 
experimentally obtained values. 


Introduction 

NASA Glenn Research Center has an ongoing research program to investigate the feasibility of developing 
jet engine fan containment systems composed of polymer matrix composite materials. To design a composite 
containment system, the ability to correctly predict the deformation and failure behavior of the composite 
under high rate loading conditions is required. Furthermore, composites composed of relatively ductile matrix 
materials are more likely to be used for containment applications. As a result, an analytical model must have 
the capability to account for nonlinearities and rate dependence in the material response. 

Previous work conducted by researchers such as Daniel, et al. (ref. 1) have determined that the material 
properties of a graphite fiber reinforced polymer matrix composite are strain rate dependent in matrix domi- 
nated deformation modes, and strain rate independent in fiber dominated deformation modes. These results 
indicate that the matrix drives any rate dependence of the composite. An example of efforts to model the rate 
dependent, inelastic response of polymer matrix composites is given by Weeks and Sun (ref. 2). In this work, 
the deformation response of the composite was modeled on the macroscopic level using a fiber orientation 
dependent plastic potential function. 

This paper describes strain rate dependent, inelastic constitutive equations for a ductile polymer and the 
incorporation of these equations into a mechanics of materials based micromechanics approach. The imple- 
mentation of the Hashin failure criterion into the micromechanics model is discussed, and efforts to implement 
the deformation model into the LS-DYNA transient dynamic finite element code are described. Verification 
studies for the deformation and ply strength models are presented for a representative composite system com- 
posed of AS-4 fibers in a PEEK thermoplastic matrix. 


Polymer Matrix Constitutive Model 

The Ramaswamy-Stouffer viscoplastic state variable model (ref. 3), which was originally developed for 
metals, was modified to simulate the rate dependent inelastic deformation of a ductile polymer. As discussed 
in (refs. 4 and 5), there are sufficient similarities between the deformation response of metals and the deforma- 
tion response of ductile polymers to permit the use of equations developed for metals to analyze ductile 
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crystalline polymers. An important point to note is that small strain conditions are currently assumed in order 
to limit the complexity of the formulation. 

In the modified Ramaswamy-Stouffer model, the inelastic strain rate, is defined as a function of the 
deviatoric stress, Sy, and back stress state variable in the form: 


e!j = D o ex P| 


] 

z 2 1 

n 

2 

[3K 2 J 



_ 


a/k7 


0 ) 


where D 0 , Z 0 , and n are material constants and K 2 is defined as follows: 


K2=^(Sij-n ij )(s ij -i2 ij ) 


( 2 ) 


The elastic components of strain are added to the inelastic strain to obtain the total strain. The following rela- 
tion defines the back stress variable rate: 




( 3 ) 


where q is a material constant, C2 m is a material constant that represents the maximum value of the back 

stress, and £<5 is the effective inelastic strain. The material constants are determined in the manner discussed 
in references 3 and 5. 

The hydrostatic stress state has been found to have a significant effect on the yield behavior of a polymer 
(refs. 4 and 6). Bordonaro (ref. 7) indicated that the proper way to account for such effects in a state variable 
constitutive model is to modify the effective stress terms. In this work, pressure dependence was included by 
multiplying the shear terms in the K 2 invariant in equation (2) by the following correction factor: 


a = 



(4) 


In this term, G m is the mean stress, J 2 is the second invariant of the deviatoric stress tensor, and (3 is a rate 
independent material constant. Since only uniaxial data was available for the polymer under consideration in 
this study, the value of the parameter (3 was determined empirically by fitting data from uniaxial composites 
with shear dominated fiber orientation angles, such as [15°]. 

To correlate the equations, the PEEK thermoplastic matrix was characterized and modeled based on 
uniaxial tensile data found in reference 7 for strain rates ranging from 1x10 6 /sec to lxlO'Vsec. While the ul- 
timate goal of this research is to model high strain-rate (up to several hundred per second) impact problems, at 
the current time high strain-rate data is not available, so all of the analyses presented in this study will be at 
relatively low strain rates. The analysis methods described in this report can be extrapolated to high strain-rate 
applications once appropriate data becomes available. The elastic modulus of the material is 4000 MPa and 
the Poisson’s ratio is 0.40. The inelastic material constants were determined to be as follows (refs. 5 and 6): 

D t) = !xl0 4 /sec, n = 0.70, Z 0 = 630 MPa, q = 310, Q m = 52 MPa, |3 = 0.45. Experimental (ref. 7) and computed 
results for strain rates of lxlO' 6 /sec and lxlO'Vsec are shown in figure 1. As can be seen in the figure, there is 
an excellent correlation between the experimental and predicted results. Note that while the tensile data was 
used to obtain the material constants, the constants were not explicitly modified in order to improve the fit of 
the computed results to the experimental data. 


Composite Micromechanics Model 

To compute the effective properties of the composite material considered in this study, a micromechanics 
technique was utilized. In micromechanics, the effective properties and response of a composite material are 
computed based on the properties of the individual constituents. The common procedure is to analyze a unit 
cell of the composite, the smallest material unit for which the response can be considered as representative of 


NAS A/TM— 1 999-20943 3 


2 



the response of the composite. For this study, the unit cell will he defined as a single fiber and its surrounding 
matrix. Furthermore, only unidirectional laminated composites at various fiber orientation angles will be 
examined. 

A mechanics of materials based micromechanics method was utilized, based on an approach developed 
by Sun and Chen (ref. 8). In this approach, the composite unit cell is divided into four subregions (fig. 2). 
Assumptions of uniform stress or uniform strain are then made within the unit cell. Along the Fiber direction 
(1 direction in fig. 2), strains in each subregion are assumed to be uniform, and the subregion stresses are 
combined using volume and local stiffness averaging. The transverse normal stresses (2 direction in fig. 2) and 
the in-plane shear stresses (1-2 direction) in subregions Af and Am are equal, while the equivalent strains are 
combined using volume averaging. The same is true for subregions B1 and B2. The total transverse (2 direc- 
tion) and in-plane shear (1-2 direction) strains in Row 1 and Row 2 (fig. 2) are assumed to be equal, and the 
equivalent stresses are combined using volume averaging. These assumptions, when combined with the consti- 
tutive equations for the Fiber and matrix, can then be used to set up systems of equations that can be solved for 
the stresses in each subregion. The subregion stresses are then used in the polymer constitutive equations to 
compute the inelastic strain rates in the matrix subregions. More details and the full equations can be found in 
reference 9. 

To verify the micromechanics equations, a composite composed of AS-4 carbon Fibers embedded in a 
PEEK matrix was analyzed. The Fiber volume fraction is 0.62, the longitudinal elastic modulus of the fiber is 
214 GPa, the transverse and in-plane shear moduli of the fiber are 14 GPa, and the longitudinal Poisson’s ratio 
is 0.2. The matrix properties are as deFmed above. In figure 3, the experimental (ref. 2) and predicted results for 
a [30°] unidirectional laminate are shown for strain rates of lxlOVsec and 0.1 /sec. As can be seen in the fig- 
ure, for both strain rates the micromechanics model predicts the rate dependence and nonlinearity of the com- 
posite response reasonably well. Note that there is some discrepancy between the experimental and computed 
results in the elastic range. As this discrepancy was not seen for [0°] or [90°] laminates, these results indicate 
that the shear moduli or transverse Poisson’s ratio used for the Fiber may not be correct. 


Ply Strength Predictions 

A variety of failure criteria exist to predict the ply level ultimate strength in polymer matrix composites. 
Several “classic” criteria are in use, which are often discussed in composite texts (e.g. (ref. 10)). Simple crite- 
ria such as Maximum Stress or Maximum Strain simply compare macroscopic (ply level) stresses (or strains) 
in each coordinate direction to maximum values. More sophisticated criteria, such as Tsai-Hill and Tsai-Wu, 
utilize quadratic (and tensor in the case of Tsai-Wu) combinations of the stresses to account for stress interac- 
tion. However, none of these criteria account for specific local failure mechanisms. 

Hashin (ref. 11) developed failure criteria that utilize quadratic combinations of the macroscopic stresses 
to approximate local failure mechanisms, such as fiber failure or matrix cracking. One advantage of utilizing 
this type of failure criteria is that by identifying specific local failure mechanisms, eventually property degra- 
dation models can be developed which allow for reductions in specific material properties and stresses as load- 
ing takes place. For this work, the Hashin failure criteria were utilized. Even though constituent level stresses 
are computed in the micromechanics, ply level strength data were more readily available and considered to be 
more reliable than constituent level strength data. Since only in-plane loads were considered in this study and 
the out-of-plane stresses were reasonably small, the plane stress versions of the criteria were employed. How- 
ever, in the future study of impact problems, where out-of-plane stresses will be signiFicant, the full three- 
dimensional version of the failure criteria will be used. Since only tensile loads were significant, only the ten- 
sile fiber failure and tensile matrix failure criteria will be presented here. For the purposes of this study, once 
either of the failure criteria was violated, property degradation was neglected, and total composite failure was 
considered to have occurred. 

Failure criteria for each of the failure modes are as follows. In each of the expressions, a ]} is the stress 
component, X T is the ply tensile strength in the longitudinal (fiber) direction, Y T is the tensile strength in the 
transverse direction, and X s is the ply shear strength. Failure is considered to occur when the value of the 
expression becomes greater than or equal to one (1). The stresses and strengths are assumed to be in the local 
material axis system, so the stresses must be transformed from the structural axis system to the material axis 
system for off-axis composites. Tensile fiber failure is predicted by using the following expression: 
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Tensile matrix failure is predicted using the following expression: 
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For the AS4/PEEK material, the longitudinal tensile strength was set to 2070 MPa, and the transverse ten- 
sile strength was set to 83 MPa (ref. 12). These values were determined to be strain rate independent for this 
material. The longitudinal tensile strength of carbon fiber reinforced composites has been found to be rate 
independent by a variety of researchers (e.g. (ref. 1)). For this material, the transverse modulus was found to be 
rate independent (ref. 2), leading to the assumption that the transverse tensile strength was most likely also 
rate independent. The in-plane shear strength was considered to be rate dependent, and was correlated based 
on data from [15°] off-axis laminates found in (ref. 2). This orientation was chosen since it is a shear domi- 
nated fiber orientation angle. From this data, the shear strength for a strain rate of lxlO’Vsec was determined 
to be 63 MPa, and the shear strength for a strain rate of 0.1 /sec was determined to be 88 MPa. The strength 
values and the experimental values for the failure stresses shown below likely have some scatter, but the sta- 
tistical data was not available. 

Using the given strength values, failure stresses were predicted for [30°] and [45°] laminates for both strain 
rates. The predicted and experimental (ref. 2) results for a strain rate of 1x10" Vsec are shown in table I, and 
the results for a strain rate of 0.1/sec are shown in table II. In all cases, failure was predicted to be due to ten- 
sile matrix failure. 

For both strain rates and fiber orientations considered, the comparison between the predicted and experi- 
mental values is quite good. In particular, the predicted results most likely fall within the scatter of the 
experimental data. The results indicate that the failure criteria are able to predict ply failure for a variety of 
fiber, orientations and strain rates. The results for this material also indicate that even when some approxima- 
tions are required in determining the ply failure stresses, reasonable results can still be obtained. 


Finite Element Implementation 

In order to simulate the impact response of a composite structure, finite element methods are required. 

With that goal in mind, the matrix constitutive equations and micromechanics model described above have 
been implemented into LS-DYNA (ref. 13), a commercially available transient dynamic finite element 
code. The implementation was accomplished through the use of the user defined subroutine option. LS-DYNA 
uses explicit central difference integration methods to integrate the rate equations. In the user defined material 
subroutine, strain increments are passed into the routine, and the total stresses at the end of the time increment 
must be computed. 

Rocca and Sherwood (ref. 14) implemented a version of the Ramaswamy-Stouffer state variable equations 
into LS-DYNA for analysis of polycarbonate. Following their formulation, equations (1) to (4) were con- 
verted into an incremental form. The back stress rate equation (eq. (3)) was modified as follows: 

AQjj = | q£2 m AeJ - qti.jAe' (7) 

where the variables are as defined earlier and the A symbol before a term indicates an increment of the vari- 
able. The back stress term fly is the total value of the back stress at the end of the previous time increment. 

The increment in inelastic strain Ae^ 1 is computed for the current time increment by taking the inelastic strain 
rate computed in the previous time increment (eq. (1)) and multiplying it by the value of the current time 
increment, which is a forward Euler type of approximation (ref. 15). The new value of the back stress for the 
current time increment is computed by adding the back stress increment to the previous value of the back 
stress. The stress and deviatoric stress increments are computed using the strain increments and the inelastic 
strain increments. The stress increments are then added to the total stresses from the previous time step to 
determine the new values of the total stresses. The inelastic strain rate for the current time step is then com- 
puted using equation (1). 
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To test the implementation of the matrix constitutive equations, the tensile response of the PEEK material 
was computed using a finite element analysis at a strain rate of 1.0/sec. The results were compared to those 
computed using the stand-alone computer code that was used to generate the results shown in figure 1 . Note, 
however, that the stresses at a strain rate of 1.0/sec are higher than those shown in figure 1. Four noded shell 
elements were used in a square mesh, ten elements on a side. One side of the model was clamped, and the 
other side had an enforced displacement to simulate a strain controlled load application. The finite element 
model was designed to simulate the behavior of the polymer at an infinitesimal material point. A mesh con- 
vergence study showed that this mesh was adequate. 

The analytical and finite element results are shown in figure 4. As can be seen in the figure, the finite 
element results compare reasonably well to those computed using the stand alone computer code. The discrep- 
ancy between the finite element and analytical results is most likely due to the nature of the integrator used 
on the rate equations. A fourth order Runge-Kutta integrator (ref, 15) was used in the stand alone computer 
code used to compute the analytical results. As mentioned earlier, a central difference explicit integrator is 
used within LS-DYNA. Furthermore, forward Euler types of approximations were used in converting the rate 
equations into an incremental form. The forward Euler method is less accurate than the Runge-Kutta method 
(ref. 15). Since the matrix constitutive equations are stiff differential equations (ref. 3), the nature of the inte- 
grator used can significantly affect the results. Future efforts may concentrate on improving the nature of the 
integration used within the user defined material subroutine. Currently, efforts are also underway to combine 
the matrix constitutive equations with the micromechanics method described above to be able to compute the 
rate-dependent nonlinear response of polymer matrix composites using the LS-DYNA finite element code. 


Conclusions 

In this study, rate-dependent inelastic constitutive equations have been developed to simulate the deforma- 
tion response of ductile polymers. The constitutive equations were implemented within a mechanics of materi- 
als based micromechanics technique. The Hashin failure criteria were utilized to predict the ply ultimate 
strengths. Preliminary efforts to implement the deformation model into the LS-DYNA finite element code were 
also discussed. Predictions made using the analytical model for an AS4/PEEK composite system compared 
well to experimental values. 

Future efforts will include completing the implementation of the deformation model into the LS-DYNA 
finite element code. Furthermore, full implementation of failure and strength models (including property deg- 
radation) into LS-DYNA will be completed. High strain rate experiments will be conducted on a representative 
polymer matrix composite, and the deformation model will be characterized and validated for high strain rate 
conditions. Furthermore, the deformation model will be extended into the large deformation regime, and the 
developed techniques will be extended to the analysis of woven composites. 
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TABLE L— FAILURE STRESS PREDICTIONS 
FOR AS4/PEEK: STRAIN RATE = lx1Q- 5 /sec 



Predicted 
failure stress, 
MPa 

Experimental 
failure stress, 
MPa 

[30°] Laminate 

130 

140 

[45°] Laminate 

98 

104 


TABLE II —FAILURE STRESS PREDICTIONS 
FOR AS 4/PEEK: STRAIN RATE = 0.1 /sec 



Predicted 
failure stress, 
MPa 

Experimental 
failure stress, 
MPa 

[30°] Laminate 

165 

170 i 

[45°] Laminate 

114 

112 
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Figure 1 : Model Correlations of PEEK Thermoplastic at Strain Rates of lxl0 -6 /sec (IE-06) and 
IxltHVsec (IE-03). 
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Figure 2: Geometry and Layout of Unit Cell Model. 
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Figure 3: Model Predictions for AS4/PEEK [30°] Laminate at Strain Rates of 1x10" ^/sec (IE-05) 
and 0. 1 sec. 



Strain 

Figure 4: Comparison of Predictions by Stand-Alone Computer Code (Analytical) and LS-DYNA 
Finite Element Analysis (FEM) for PEEK Thermoplastic at Strain Rate of 1.0/sec. 
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